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-^A  relative  waveform  inversion  method  has  been  developed  which 
determines  models  that  correctly  predict  waveform  changes  in  a 
set  of  observed  seismograms.  By  modeling  changes  in  waveforms 
rather  than  absolute  waveforms,  a  number  of  relatively  poorly 
determined  quantities,  such  as  receiver  function  and  Q,  may  be 
eliminated  from  the  problem  and  greater  resolution  of  key  model 
parameters  obtained.  The  method  has  been  applied  to  the  estimation 
of  overall  amplitude.  pP  amplitude  and  pP  delay  time  for  a  suite  of 
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of  azimuthal  variation  in  pP.  The  existence  of  events  with  radical 
different  pP  amplitudes  suggests  that  biases  may  exist  in  current 
relative  yield  estimation  techniques,  an  effect  which  is  corrected 
for  in  the  relative  waveform  method 
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I .  INTRODUCTION 

Two  basic  classes  of  methods  have  generally  been  used 
in  the  analysis  of  seismic  data.  The  older,  and  still  most 
prevalently  used,  class  of  methods  we  will  call  data  para- 
metrization  methods.  The  other,  more  recently  developed 
class  of  methods,  we  will  call  modeling  methods. 

Data  parametrization  methods  are  those  techniques  in 
which  we  attempt  to  find  models  that  fit  only  selected 
aspects  of  the  data.  That  is,  a  parametrization  of  the  data 
is  chosen  which  does  not  completely  describe  the  data,  and 
these  parameters  are  then  used  to  determine  suitable  models. 

Examples  of  such  par^metrizations  include  travel  times, 
spectral  slopes  and  power  spectra.  While  these  features  of 
the  data  may  contain  much  of  the  information  concerning  the 
model  parameter  of  interest,  they  are  not,  in  themselves, 
sufficient  to  completely  recover  the  data. 

In  modeling  studies,  on  the  other  hand,  no  parametri¬ 
zation  of  the  data  is  used,  although  data  may  be  prefiltered 
and  windowed.  Instead,  we  attempt  to  construct  best  fitting 
estimates  of  the  prefiltered  data,  using  estimates  of  all 
necessary  model  parameters. 

The  major  advantage  of  data  parametrization  methods  is 
that  it  is  often  possible  to  find  simple  functional  forms 
relating  the  parametrized  data  to  the  parametrized  model. 

This  allows  the  inverse  problem  to  be  done  quite  simply, 
often  as  a  simple  linear,  least  squares  problem  even  when  no 
such  simple  solution  of  the  full  inverse  problem  exists. 
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A  second  advantage  is  the  ease  with  which  relative  varia¬ 
tions  in  a  single  model  parameter  may  be  studied.  A  variety 
of  spectral  ratio  and  similar  techniques  are  available  which 
allow  a  single  or  small  number  of  model  parameters  to  be 
isolated  and  variations  in  these  data  characteristics  to 
completely  define  the  model  studied  independently  of  other 
factors  necessary. 

The  principal  advantage  of  modeling  studies,  on  the 
other  hand,  is  the  greater  resolution  that  can  often  be 
obtained  by  making  use  of  non-parametrized  data  seismograms. 
One  example  of  this  approach  is  the  fine  structure  obtained 
through  waveform  inversion  of  triplication  data,  compared  to 
the  models  obtainable  through  more  standard  travel-time 
analysis  (Mellman,  1980). 

Use  of  the  additional  data  in  the  non-parametrized 
seismogram  means  that  it  is  always  possible  to  find  a  model¬ 
ing  method  that  gives  greater  model  resolution,  or  smaller 
model  variance,  than  is  obtainable  in  a  corresponding  data 
parametrization  study.  Unfortunately,  it  is  often  not 
possible  to  achieve  that  resolution  in  practice.  Such 
results  can  only  be  achieved  through  the  proper  choice  of 
error  function,  a  fact  which  is  ignored  in  most  model  stud¬ 
ies.  Indeed,  most  forward  modeling  studies  rely  on  fitting 
by  eye,  rather  than  on  optimization  of  an  error  function. 
While  this  gives  the  researcher  considerable  flexibility  in 
choice  of  which  aspects  of  the  data  are  important,  it  is 
difficult  to  assess  how  well  the  model  is  constrained. 
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During  the  course  of  a  study,  one  usually  achieves  a  strong 
intuitive  feel  for  model  constraint,  but  it  is  virtually 
impossible  to  pass  this  intuition  on  to  others  in  the  field. 

Moreover,  due  to  the  existence  of  unrecognized  tradeoffs, 
this  intuitive  feel  for  the  model  constraints  may  be  mis¬ 
leading  or  incorrect. 

The  use  of  formal  inverse  methods  in  modeling  studies 
requires  the  adoption  of  an  objective  error  function.  Even 
in  this  approach,  however,  little  attention  is  paid  to 
choosing  error  functions  that  provide  the  best  constrained 
models.  Proper  choice  of  these  error  functions  is  in¬ 
timately  related  to  effects  of  various  types  of  noise  on  the  ! 

final  solution.  This  includes  both  modeling  noise,  which  is 
almost  universally  ignored,  and  the  more  commonly  considered 
data  noise.  This  topic  is  relatively  poorly  understood,  and 
much  additional  work  is  necessary,  paritcularly  for  non¬ 
linear  problems. 

Even  without  optimal  choice  of  an  error  function, 
modeling  studies  can  often  indicate  the  presence  of  noise 
problems  in  cases  where  data  parametrization  methods  give  no 
indication  of  any  problems.  This  is  especially  true  in  the 
case  of  modeling  noise.  By  modeling  noise,  we  mean  noise 
induced  by  incomplete  or  incorrect  model  parametrization. 

It  is  often  only  by  detailed  analysis  of  data  misfits  in 

modeling  studies  that  it  is  possible  to  decide  that  a  more 

general  model  is  appropriate.  This  can  effect  not  only 

estimates  of  the  model,  but  estimates  of  the  uncertainty  of 

model  parameters  as  well.  I 
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A  major  disadvantage  of  modeling  studies  has  been  the 
absence  of  a  stable,  generally  applicable  modeling  method 
equivalent  to  spectral  ratio  methods  for  the  study  of  rel¬ 
ative  changes  of  a  suite  of  seismograms.  Such  methods  do 
exist  for  specific  cases,  such  as  attenuation  studies 
(Burdick  and  Helmberger,  1974)  where  group  properties  of  the 
attenuation  operator  may  be  exploited  to  provide  stable 
estimates  of  differential  attenuation.  In  general,  however, 
the  only  method  available  for  relative  studies  has  been 
deconvolution,  which  is  notorious  for  its  instability  under 
even  the  best  of  circumstances. 

What  criteria  would  we  like  in  a  relative  modeling 
method?  First,  we  require  the  method  to  be  symmetric  with 
respect  to  the  data.  That  is,  the  results  should  be  in- 
pendent  of  the  designation  of  any  particular  datum  as  the 
reference  datum.  Thus  no  datum  is  given  a  preferred  pos¬ 
ition  in  the  method.  Second,  the  method  should  be  stable. 
Third,  we  would  like  to  have  sufficient  local  linearity  to 
make  an  iterative  linearized  inversion  procedure  feasible. 
We  note  that  deconvolution  methods  satisfy  none  of  these 
criteria.  In  the  following  section,  we  describe  a  method 
that  does  satisfy  these  criteria  and  its  inverse  as  applied 
to  seismic  waveform  data.  We  will  call  this  method  relative 
waveform  inversion. 
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II.  RELATIVE  WAVEFORM  INVERSION 

Consider  two  seismograms,  f1(t)  and  f2(t)  composed  of  a 
common  transfer  function  T(t)  and  a  portion  S^(t)  which 
differs  from  seismogram  to  seismogram.  If  we  choose  f^(t) 
and  f2(t)  t0  seismograms  rcorded  at  a  single  station  in 
the  far  field  for  two  events  in  a  common  source  region  then 
the  transfer  function  T  would  contain  instrument,  path,  and 
surce  properties,  and  receiver  response  functions  common  to 
both  seismograms.  The  remaining  source  and  propagation 
terms  are  included  in  S^.  Then,  by  definition 


f^t)  =  T(t)  *  Si(t) 


If  this  expression  is  exact,  then 


f.(t)  *  S,(t)  =  f„(t)  *  S. (t) 


This  suggests  that  minimization  of  the  error  functional 


e12  "  I  I  W 

F12(t)  =  £1*s2-f2*s1 


where  ( 3 ) 


I  X(t)  | 


|  X(t)  *  W(t)  |  2dt. 


will  provide  a  good  estimate  of  the  functions  S^. 

The  weighting  function  W  acts  as  a  prewhitening  filter 
and  is  chosed  to  extend  the  usable  signal  bandwidth  while 
still  rejecting  noise.  In  general,  choice  of  W  will  depend 
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on  the  frequency  content  of  f^  and  and  on  the  noise 
spectra.  For  the  applications  discussed  below,  simple  fun¬ 
ctions  of  the  form  W(w)  =  (l+kiu)  were  found  to  be  adequate. 
Choice  of  the  error  function  F12(t)  in  equation  3  leads  to  a 
relative  waveform  method  that  satisfies  the  three  criteria 
presented  in  the  introduction.  The  error  functional  e12  is 
invarient  to  interchange  of  indicies,  assuring  that  neither 
seismogram  occupies  a  preferred  position.  Further,  the 
error  function  F12  consists  of  differences  of  the  con¬ 
volutions  of  well-behaved  functions  which  assures  stability. 
While  the  local  linearity  of  F12  with  respect  to  model  pa¬ 
rameters  depends  on  the  choice  of  model  parametrization,  the 
functional  form  of  F12  gives  promise  of  sufficient  linearity 
for  a  wide  variety  of  model  parameterization  to  make  iter¬ 
ative  linearized  inversion  feasable. 

The  error  functional  e12  in  equation  2  may  be  readily  ex¬ 
tended  to  multiple  recording  stations.  Let  e112  be  the 
error  functional  for  the  i1*1  station  with  fj^  the  seis¬ 
mogram  for  event  j  at  station  i.  We  may  then  define  a  new 
error  functional . 

<ei2>2  = 1  vi<eL>2  (4) 

with  vi  =  (/f^tjdtr*5  ( /^.(tjdtr4* 

Here  is  a  station  weighting  factor,  used  to  prevent  high 
amplitude  observations  from  dominating  the  solution. 
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One  obvious  application  of  the  relative  waveform  in¬ 
version  technique  is  the  estimation  of  "depth  corrections" 
to  the  yields  of  underground  nuclear  explosions.  It  has 
been  appreciated  for  some  time  that  surface  reflections  have 
a  major  effect  on  both  waveform  and  amplitude  of  shallow 
events.  Efforts  to  account  for  thses  effects  have  often 
taken  form  of  modeling  studies,  where  elastic  wave  veloc¬ 
ities  are  used  to  predict  arrival  times  and  amplitudes  of 
surface  reflections.  See,  for  example,  Burdick  and  Helm- 
berger  (1979).  Such  approaches  have  been  criticized  for  the 
linear  assumptions  used  in  what  is  clearly  a  nonlinear 
regime.  In  particular,  it  is  often  observed  that  arrival 
times  of  surface  reflections  are  delayed  with  respect  to 
predictions  made  on  the  basis  of  know  elastic  velocities. 
This  is  true  both  observational ly  and  theoretically  based  on 
non-linear  finite  difference  calculations  (Trulio  et  al, 
1980;  Mellman  et  al,  1980).  In  addition,  several  workers 
(Bache  et  al,  1979;  Blandford  et  al,  1979)  have  presented 
observational  evidence  that  at  least  at  high  frequencies, 
amplitudes  of  surface  reflections  can  be  greatly  reduced 
compared  to  estimates  based  on  linear  elastic  theory.  These 
amplitude  effects  can  be  seen,  once  again,  in  theoretical 
non-linear  models.  We  would  like  any  inversion  technique  to 
include  these  effects. 

While  anomalous  amplitude  and  travel  times  of  surface 
reflections  can,  to  some  extent,  be  estimated  using  tradi¬ 
tional  modeling  techniques  (Burdick  and  Helmberger,  1979) 
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uncertainties  in  time  function  and  attenuation  operators 
limit  the  resolution  obtainable.  In  such  circumstances,  we 
expect  that  relative  waveform  methods  will  be  of  maximum 
utility. 

If  we  restrict  ourselves  to  consideration  of  events  of 
approximately  the  same  size  in  the  same  source  region  (and 
thus,  presumably  in  the  same  geologic  material),  recorded  at 
the  same  teleseismic  station,  it  is  reasonable  to  assume 
that  the  major  differences  in  seismograms  for  different 
events  will  be  caused  by  an  overall  scale  factor  change  and 
changes  in  pP  travel  time  and  amplitude.  We  therefore 
choose  the  common  portion  of  the  synthetic  seismograms,  T(t) 
in  equation  2,  to  include  the  normalized  source  time  func¬ 
tion,  geometric  spreading,  anelastice  attenuation,  receiver 
function  and  instrument  transfer  function.  The  functions 
S^(t)  are  of  the  form 

51  =  cx  (6(t)  -  ai  6(t-ti))  (5) 

52  =  c2  (6 (t-t0  )-a2  6  ( t-t0-t2 ))  . 

In  order  to  reject  the  trivial  minimum  of  error  functional 
e12  in  equation  4,  namely  cx  =  c2  =  0,  we  impose  the  con¬ 
straint 


Ci  •  c2  =  1  or  equivalently  (6) 

Ci  =  1/c 
C2  ”  c. 

For  parametrization  of  more  than  two  S^,  we  procede  as 
above,  with  the  constraint  Ilc^  =  1.  In  this  case,  we  min¬ 
imize  a  new  error  function 
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i>j  c  c 

i  j 


(7) 


For  simplicity,  however,  we  will  restrict  the  following 
discussion  to  the  two  event  case. 

At  this  point,  it  is  convenient  to  define  a  new  normal¬ 
ized  concatinated  function  En.  This  is  formed  by  con- 
catinating  functions  V.  ^F1i2*W  where  W' (u))=l+kw.  The  con¬ 
catinated  functions  are  then  sampled  at  discrete  time  points 
.  Thus  En  is  the  n^  time  point  of  the  concatinated  pre¬ 
filtered  and  normalized  error  functions.  We  note  that 
minimization  of 

e2  h  I(En)2  (8) 

is  equivalent  to  minimization  of  the  error  function  e12,  but 
leads  to  a  somewhat  more  simple  formulation. 

The  actual  determination  of  the  is  carried  out  using  an 
iterative  stabilized  linearized  least-squares  inverse  me 
thod.  We  let  nij  be  the  model  parameters  that  define  all  the 
s^.  Then  some  change  in  model  parameters  6mj  will  cause 
some  change  6En  in  the  error  function  ER.  We  wish  to  find 
6m j  to  satisfy  Ei  +  6Ei  =  0  in  a  least  squares  sense.  By 
Taylor  expansion,  we  find 
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6Ei  =  ffi  fimj  =  j  6nij .  (9) 

anij 

This  gives  us  the  linearized  equation 

Aij  6mj  +  Ei  =  0.  (10) 


The  weighted  least  squares  solution  to  this  equation  is 
given,  in  matrix  notation,  by 


6m  =  -(ATA  +  yD"1  ATE 


(11) 


T 

where  y  -  a  trace  (A  A)  xs  a  damping  factor  used  to  sta¬ 
bilize  the  inversion.  The  form  of  A,  the  derivative  matrix, 

is  given  in  the  appendix.  The  inversion  procedure  is  ap- 

k+1 

plied  iteratively,  so  that  the  model  estimate  m^  after 
k+1  iterations  is  given  by 


mjktl  =  mjk+1  *  8mjk+l' 


(12) 


We  note  that,  due  to  the  non-linearity  of  the  problem,  the 
matrix  A^j  must  be  recomputed  at  each  iteration. 

The  form  of  the  error  function  in  equation  3  and  the 
model  parametrization  used  give  the  relatvie  waveform  in¬ 
version  several  interesting  properties.  In  general,  dif¬ 
ferences  in  pP  times  and  amplitudes  for  two  events  will  be 
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much  better  constrained  than  the  amplitude  or  travel  time  of 
and  individual  pP.  As  an  example,  consider  the  case  in 
which  seismograms  for  two  events  are  identical.  Thus 

fi  (t)  =  f2(t)  and 

S^ * f 2 “ S2 * f x =0  if  Sx=S2*  (13) 

In  this  case,  there  are  no  constraints  on  individual  pP 
times  or  amplitudes  at  all,  while  time  and  amplitude  dif¬ 
ferences  for  the  two  events  remain  well-constrained.  Reso¬ 
lution  of  individual  times  and  amplitude  improves  as  two 
events  become  increasingly  dissimilar.  This  fact  is  borne 
out  by  examination  of  the  model  varience  matricies  for 
synthetic  and  actual  data  runs. 

Another  interesting  aspect  of  relative  waveform  in¬ 
version  is  the  ability  to  determine  pP  arrival  times  even 
when  these  times  are  less  than  the  inverse  of  the  high 
frequency  cutoff  for  bandpassed  data.  This  is  due  to  the 
fact  that  the  norm  of  the  error  function,  viewed  in  the 
frequency  domain,  is  minimized  by  matching  phase  at  all 
available  frequencies.  This  is  in  contrast  to  cepstal 
methods,  for  instance,  which  require  that  two  arrivals  have 
a  n  phase  shift  at  some  frequency  in  order  to  produce  an 
arrival  time  estimate. 

If  this  frequency  is  high  enough  to  be  outside  usable 
signal  passband,  or  if  the  second  arrival  is  not  a  spectral 
shadow  of  the  first  arrival  at  high  frequencies,  the  cep- 
stral  method  will  fail.  Under  these  circumstances,  the 
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relative  waveform  inversion  will  still  produce  a  reliable 
estimate  valid  in  the  data  passband. 

Although  much  of  the  development  in  this  section  has 
been  done  in  terms  of  a  specific  problem,  the  basic  relative 
waveform  method  has  far  more  general  applicability.  The 
error  function  defined  in  equation  3  and  its  rather 
straightforward  generalization  to  multiple  receivers  and 
many  events,  is  in  no  way  tied  to  the  specific  source 
estimation  problem  in  question.  In  particular,  one  example 
in  which  the  relative  waveform  method  is  likely  to  yield 
good  results  is  in  the  estimation  of  earth  structure 
controlling  triplications.  Conventional  modeling  methods  are 
quite  powerful  in  this  type  of  problem,  but  suffer  from  the 
necessity  of  obtaining  good  source  and  attenuation 
estimates,  even  though  both  source  and  attenuation  are 
stationary  over  the  distance  range  of  interest.  Thus,  such 
studies  often  must  be  preceded  by  detailed  source  studies, 
while  use  of  short  period  waveform  data  is  precluded  due  to 
the  lack  of  adequate  source  models  at  short  periods.  By 
using  the  relative  waveform  method  for  all  possible  pairs  of 
observations  for  a  given  event,  and  including  a  number  of 
events,  move  out  and  relative  amplitude  information  may  be 
recovered  from  waveforms  without  the  necessity  of  detailed 
source  studies. 

As  in  any  other  method,  the  validity  of  the  results  of 
a  relative  waveform  inversion  depends  on  the  validity  of  the 
assumptions  that  go  into  the  method.  We  would  therefore 
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like  to  delineate  the  assumptions  made,  and  to  explore  the 
effects  of  these  assumptions,  where  possible. 

For  tightly  grouped  events  observed  at  a  single  sta¬ 
tion,  it  is  reasonable  to  assume  that  attenuation,  instru¬ 
ment  transfer  function,  and  receiver  function  will  be  common 
for  all  events.  Use  of  a  single  source  time  function  for  all 
events,  however,  is  less  easily  justified.  It  is  certainly 
unjustified  for  events  of  greatly  different  size  or  events 
occuring  in  grossly  different  materials,  such  as  granite  and 
unsaturated  sediments.  If,  however,  the  events  are  in 
similar  materials,  it  is  anticipated  that  the  degree  of 
overshoot  in  the  time  functions  of  the  two  events  will  be 
similar  (Haskell  1967).  Also,  for  small  to  intermediate 
events,  changes  in  the  rise  time  of  the  source  time  function 
will  produce  only  minor  effects  in  the  final  seismogram. 
This  has,  to  some  degree,  been  tested  using  synthetic  source 
functions  derived  from  finite  difference  calculations  by 
Perl  and  Trulio,  1981,  for  a  variety  of  events.  Results  of 
synthetic  studies  on  source  function  bias  in  relative  wave¬ 
form  inversions  using  von  Seggern-Blanford  sources  are 
presented  in  a  later  section, 

A  second  cause  for  concern  lies  in  the  use  of  a  single, 
frequency  independent  arrival  to  represent  all  of  the  non¬ 
constant  effects  of  near  source  propogation.  Certainly,  in 
the  presence  of  complex  crustal  structure,  small  changes  in 
source  location  can  result  in  differing  numbers,  amplitudes 
and  arrival  times  for  crustal  reflections.  Also,  as  men- 
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tioned  earlier,  evidence  exists  for  the  frequency  dependence 
of  the  amplitude  of  the  surface  reflection. 

Both  time  function  and  propagation  effects  can  easily 
be  explicitly  included  in  the  inversion  model.  Un¬ 
fortunately,  the  limited  bandwidth  available  in  observed 
seismograms  makes  it  doubtful  that  any  of  these  effects 
could  be  differentiated  from  a  simple  change  in  pP  amplitude 
or  arrival  time.  For  this  reason,  it  was  decided  to  allow 
all  variation  to  occur  in  the  surface  reflection,  and  to 
investigate  resolvability  of  additional  early  arrivals  at  a 
later  date. 

How,  then,  are  we  to  interpret  the  results  of  a  rela¬ 
tive  waveform  inversion?  Clearly  the  second  arrival  in  a 
narrow  frequency  band  a  time  corresponding  to  some  weighted 
average  of  these  arrival  times.  To  the  extent  that  these 
arrivals  are  dominated  by  pP,  the  second  arrival  time  and 
amplitude  will  be  accurate  estimates  of  pP  delay  and  average 
amplitude  within  a  narrow  spectral  window.  The  estimated 
overall  amplitude  factor,  however,  can  be  expected  to  remain 
valid  independent  of  the  exact  nature  of  the  constituents  of 
the  second  arrival.  Alternatively,  if  the  results  of  the 
relative  waveform  inversion  may  be  viewed  as  providing  a 
measure  of  the  degree  of  difference  of  observed  seismograms, 
expressed  in  terms  of  a  simple,  plausible  source  model.  At 
the  very  least,  this  provides  a  mechanism  for  identifying 
events  that  are  anomalous,  or  at  least  show  large  dif¬ 
ferences  from  other  events  in  ways  that  cannot  be  explained 
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by  very  small  changes  in  pP  time  and  amplitude.  And  again, 
the  important  fact  in  adjusting  yield  estimates  is  the  size 
and  timing  of  the  effective  second  arrival,  not  the  physical 
causes  of  this  arrival. 

An  additional  assumption  of  source  isotropy  is  made 
when  data  from  several  stations  at  widely  differing  azimuths 
is  used  in  the  same  inversion  run.  In  cases  of  large  scale 
tectonic  release,  asymmetric  spall,  or  lateral  variation  in 
structure,  the  effective  source  function  s^  for  a  given 
event  may  show  significant  azimuthal  variation.  A  com¬ 
parison  of  errors  for  an  inversion  containing  data  sampling 
a  range  of  azimuths  with  inversions  contain  data  sampling 
single  azimuths  should  provide  a  means  of  assessing,  and 
studying,  azimuthal  source  variations. 
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III.  SYNTHETIC  RESULTS 

A  series  of  tests  on  synthetic  data  were  conducted  to 
test  the  relative  waveform  inversion  methodology  and  its 
sensitivity  to  variations  in  input  parameters,  in  data  and 
source  characteristics,  and  in  inversion  parameters.  Some 
results  of  the  first  set  of  test  cases  are  illustrated  in 
Figures  1  and  2.  These  experiments  examined  the  sensitivity 
of  the  technique  to  choice  of  starting  values  (with  respect 
to  the  true  values)  and  to  selection  of  the  damping  para¬ 
meter  in  the  actual  inverse  procedure.  Two  synthetic  short 
period  P  wave  seismograms  were  constructed  for  these  tests. 
Both  seismograms  have  identical  instrument,  attenuation,  and 
source  time  function  characteristics,  differing  only  in  the 
time  delay  and  amplitude  of  the  second  (pP)  arrival  (see 
Figure  1).  In  order  to  examine  the  effect  the  choice  of 
starting  model,  the  inverse  assumed  initially  that  both 
seismograms  had  the  same  pP  delays  and  amplitudes  and  iter¬ 
ated  from  that  point.  In  addition,  this  procedure  was 
repeated  using  a  different  value  of  the  damping  parameter  a. 
The  results  of  both  tests  are  shown  in  Figure  2.  We  see  that 
in  both  of  these  particular  examples,  fairly  rapid  con¬ 
vergence  to  the  correct  values  was  attained.  Figure  3  com¬ 
pares  the  derived  waveforms  following  the  inverse  procedure. 
The  rate  of  convergence  is,  of  course,  controlled  by  the 
inversion  damping  parameter.  Although  in  the  present  exam¬ 
ple  both  a  values  (0.01  and  0.001)  permitted  stable  con¬ 
vergence,  further  testing  has  shown  that  the  choice  of  0.01 
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SYNTHETIC  SEISMOGRAMS 


Figure  1.  Two  test  synthetics  were  constructed  to  check  the 
convergence  of  relative  waveform  inversion.  The 
synthetics  have  the  same  instrument  response, 
attenuation  and  time  function  convolved  together 
with  different  pP  spike  arrivals.  T#1  has  a  dP  spike 
arrival  at  .30  sec.  with  an  amplitude  of  .75  (Rel. 
to  P)  while  T#2  has  a  pP  arrival  at  .20  sec.  with 
an  amplitude  of  .50. 
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Figure  2 
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Using  our  two  synthetic  seismograms  we  invert  for  pP 
arrival  time  and  amplitude  given  a  starting  estimate . 
In  this  graph  of  pP  amplitude  versus  arrival  time  we 
have  selected  our  starting  position  at  .6  and  .4  sec. 
The  inversion  successfully  converged  from  this  point 
to  correct  synthetic  amplitudes  and  arrival  times.  We 
ran  the  inversion  with  two  different  damping  factors 
to  illustrate  convergence  stability. 
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CONVERGENCE  PROPERTIES  OF  THE 
RELATIVE  WAVEFORM  INVERSION  TECHNIQUE 


SYNTHETIC  1 


pP  ARRIVAL  TIME  (sec) 


Figure  3: 

In  this  figure  we  have  the  inversion  results  from  our  test 
synthetics.  We  have  inverted  for  two  spike  seismograms  SI 
and  S2  which  converged  to  the  source  terms  convolved  into  our 
synthetics.  Given  one  synthetic  then,  Fl  we  can  construct  a 
duplicate  of  Fl  by  convolving  the  transfer  function  (S1/S2) 
with  F2.  Fl  and  F2*  (S1/S2)  are  shown  here  for  comparison. 

If  the  inversion  has  converged  to  the  exact  solution,  which  it 
has  in  this  case,  then  these  two  seismograms  will  be  identical. 
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is  highly  preferred.  While  this  value  results  in  somewhat 
slower  convergence,  it  is  always  a  stable  procedure.  While 
the  smaller  damping  afforded  by  a  a  =  0.001,  causes  gen¬ 
erally  faster  convergence,  frequently  the  smaller  damping 
causes  either  divergent  solutions  or  a  "hem-stitching" 
behavior  in  the  iterative  procedure. 

The  second  series  of  tests  with  synthetic  data  examined 
the  effects  of  random  noise  in  the  data  as  well  a^  testing 
for  distortion  of  the  final  results  by  such  factors  as 
complex  receiver  functions.  Figures  4  and  5  show  the  basic 
synthetic  data  employed,  the  same  data  with  10%  instrument 
bandpassed  noise  added,  and  the  synthetic  data  convolved 
with  observed  receiver  functions  for  the  SRO  stations  CHTO 
and  ANMO.  The  choice  of  instrument  bandpassed  noise  was 
felt  to  be  a  most  severe  test  due  to  the  fact  that  this 
noise  has  a  similar  frequency  content  to  the  signals.  The 
convergence  path  results  of  two  test  inversions,  utilizing 
different  starting  models,  of  the  synthetic  data  with  noise 
added  are  shown  in  Figure  6.  Both  inversions  converged  to 
the  same  final  values  which  differ  slightly  from  the  true 
parameters.  However,  these  noise-induced  errors  were  gen¬ 
erally  very  small  {less  than  10%  in  pP  amplitude  and  less 
than  0.1  second  in  pP  delay  time). 

Since  the  near-receiver  crustal  response  is  a  common 
feature  in  the  data,  theoretically  it  should  have  no  in¬ 
fluence  upon  a  relative  waveform  inversion.  However,  to 
test  this  basic  assumption,  especially  in  light  of  the 
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Figure  4. 


Shown  above  are  four  basic  synthetic  seismograms  used 
to  simulate  one  of  the  two  station  seismograms  needed 
for  relative  waveform  inversion.  The  first  synthetic 
has  only  an  instrument,  source  time  function,  attenua¬ 
tion  operator  and  spike  (P,pP)  seismogram  convolved 
together.  In  addition  the  second  trace  has  10%  random 
instrument  filtered  noise  added  to  it.  Synthetics  3  and 
4  have  know  receiver  functions  for  SRO  stations  CHTO 
and  ANMO  convolved  in  instead  of  noise. 
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Figure  5.  Shown  above  are  four  basic  synthetic  seismograms  used 
to  simulate  one  of  the  two  station  seismograms  needed 
for  relative  waveform  inversion.  The  first  synthetic 
is  composed  of  an  instrument,  source  time  function, 

♦  attenuation  operator,  and  spike  (P,pP)  seismogram. 

In  addition  to  these  convolution  terms,  the  second 
trace  has  10%  random  instrument  filtered  noise  added 
to  it.  Synthetics  3  and  4  have  known  receiver  functions 
for  SRO  stations  CHTO  and  ANMO  convolved  in  instead 
of  noise. 
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Figure  5 


Figure  6: 

Relative  waveform  inversion  was  performed  on  the  two  test 
synthetics  with  random  10%  instrument  bandpassed  noise  added 
to  each  signal.  Using  the  same  two  starting  models  we  plotted 
the  inversion  path  in  term  of  Pp  amplitude  versus  arrival  time. 
The  maximum  error  in  this  test  was  only  10%  in  amplitude  and 
0.1  seconds  which  may  be  reduced  (by  )  by  adding  more 

stations  to  the  inversion. 
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occasionally  very  complex  nature  of  the  near-receiver  ef¬ 
fects,  the  influences  of  the  receiver  effects  at  the  SRO 
stations  CHTO  and  ANMO  (Figures  4  and  5)  were  specifically 
tested.  Since  both  of  these  stations  were  used 
extensively  in  the  applications  of  this  technique  to  ob¬ 
served  data,  this  was  felt  to  represent  a  particularly 
important  test.  Figures  7  and  8  plot  the  convergence  paths 
for  the  inversion  test  using  the  ANMO  and  CHTO  receiver 
functions.  A  maximum  error  of  6%  was  observed. 

In  all  of  the  tests  conducted,  the  errors  introduced  by 
outside  factors  were  encouragingly  small  in  relative  P  wave 
amplitude,  pP  amplitude,  and  pP  delay  time.  Moreover,  even 
these  relatively  small  errors  can  be  further  reduced  by 
utilizing  more  than  two  stations  in  the  inversion  procedure. 

The  final  series  of  synthetic  data  tests  examined  the 
influences  of  source  time  function  on  the  technique.  The 
assumption  that  the  source  time  function  is  identical  amoung 
the  events  examined  would  seem  to  be  the  most  questionable, 
and  these  part  of  the  study  attempts  to  quantify  the  ef¬ 
fective  bias,  if  any,  introduced  by  variations  in  source 
time  function.  Synthetic  data  were  constructed  using  von 
Seggern  and  Blandford  (1972)  source  functions  with  a  dis¬ 
tribution  of  overshoot  and  rise  time  characteristics.  The 
von  Seggern  -  Blandford  source  time  function,  in  the 
far- field,  is  given  by: 
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Figure  7.  The  observed  receiver  function  for  SRO  station  ANMO  j 

was  convolved  with  both  our  test  synthetics  before 
inversion.  Shown  here  are  the  convergence  paths  from 
two  different  starting  models  plotted  in  pP  amplitude- 
I  arrival  time  space.  Station  ANMO  has  a  slightly  greater 

effect  on  convergence  than  did  CHTO  no  doubt  because 

it  acted  as  a  lowpass  filter  and  eliminated  necessary 

high  frequency  waveform  details.  i 

3 


AMPLITUDE 


ANMO  RECEIVER  FUNCTION 
BIASING  TEST 


Pp  ARRIVAL  TIME  (SEC) 


Figure  7 


33 


figure  8.  A  receiver  function  for  the  SRO  station  CHTO  was 
convolved  with  our  two  test  synthetics  before  in¬ 
version  to  check  common  filtering  effects.  Shown 
above  are  the  inversion  paths  from  two  different 
starting  models  plotted  in  terms  of  pP  amplitude 
versus  arrival  time.  As  expected  the  common  con¬ 
volution  factor  has  very  little  effect  on  converg¬ 
ence. 
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Figure  8: 
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V(t)  =  yq  ke"kT  (  (l+2B)kt  -  Bk2  x2  ) 

where  k  is  the  risetime  parameter,  B  controls  the  degree  of 
overshoot,  and  y0  is  an  overall  sealing  factor.  We  shall  be 
only  concerned  here  with  examining  the  effects  of  variations 
in  k  and  B. 

Typical  values  of  k  and  B  for  an  event  with  a  yield  on 
the  order  of  100  kt  are  10  and  1  respectively.  For  the  pur¬ 
poses  of  these  inversion  tests,  these  values  were  selected 
as  the  mean  and  variations  about  these  values  were  used  to 
construct  the  test  dataset.  Graphical  representations  of 
the  source  time  functions  with  k  varied  between  5  and  15  and 
B  varied  between  0  and  3  are  shown  in  Figure  9.  The  syn¬ 
thetic  siesmograms  were  constructed  by  convolving  the  time 
functions  shown  in  Figure  9  with  a  short  period  SRO  instru¬ 
ment  response,  an  attenuation  operator,  6-function  P  and  pP 
arrivals  and  CHTO  and  ANMO  receiver  functions .  The  final 
synthetics  are  shown  in  Figures  10  and  11.  In  all  cases, 
the  pP  arrival  was  inserted  with  a  .4  second  time  delay  and 
an  amplitude  of  .8  relative  to  the  direct  P  wave. 

A  series  of  two  station  inversion  tests  were  conducted 
using  the  standard  (k=10,  B=l)  synthetic  seismogram  and  the 
other  synthetic  data.  The  starting  point  in  each  test  was 
the  true  values  (pP  amplitude  =  .8;  pP  delay  =  .4  sec.).  One 
important  factor  noted  here  and  in  all  previous  tests  was 
that  the  inversion  becomes  more  constrained  and  resolution 
increases  apparent  differences  in  the  observed  waveforms 
increase . 
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Figure  9a: 

Source  time  functions  from  the  Von  Seggern  and  Blandford  model 
are  shown  above  for  a  range  of  risetimes  with  fixed  overshoot 
(B=l) .  From  these  time  functions  we  construct  station  synthetics 
and  test  for  inversion  bias  due  to  a  variable  source  term. 

Figure  9b: 

Source  time  functions  shown  above  cover  a  range  of  overshoot 
with  the  risetime  fixed  at  K=10. 
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i  Figure  10a: 

Synthetic  seismograms  are  shown  for  a  range  of  overshoot  parameter 
B  and  constant  rise  time  K.  Receiver  function  appropriate  for 
CHTO  has  been  used. 

Figure  10b: 

As  in  Figure  10a,  but  with  ANMO  receiver  function. 
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Figure  11a.  Synthetic  siesmograms  are  shown  for  a  range  of 

risetime  parameters  K  with  constant  overshoot  B. 
Receiver  function  for  CHTO  has  been  used. 


Figure  lib.  As  in  figure  11a  buth  with  ANMO  receiver  function. 


Figure  11a:  Figure  lib 
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We  first  consider  the  effects  of  variations  in  over¬ 
shoot  parameter.  Figure  12  illustrates  the  biasing  effect 
of  changes  in  overshoot  on  the  results  of  the  relative 
waveform  inversion.  In  general,  the  derived  pP  arrival  time 
(Figure  12b)  is  little  effected  by  changes  in  overshoot. 
Only  in  the  extreme  case  of  a  0:1  ratio  in  the  B  value 
between  the  two  seismograms  does  the  error  exceed  .02  sec¬ 
onds  and  even  then  the  error  reaches  only  0.07  seconds.  The 
derived  pP  amplitude  is  somewhat  more  sensitive  (Figure  12a) 
although  again  the  bias  exceeds  about  15%  only  when  the 
ratio  of  B  values  becomes  extreme  (0:1  or  3:1).  In  prac¬ 
tice,  such  changes  are  unlikely  to  be  encountered  in  the 
observed  data  sets  we  contemplate. 

Figure  13  summarizes  the  effects  of  variations  in 
risetime  (k)  on  the  method.  Here  we  find  somewhat  of  a 
reversal  in  trend  in  comparison  to  the  time  and  amplitude 
variations  produced  by  changes  in  the  B  value.  The  pP 
amplitude  is  almost  insensitive  to  variations  in  k  except  at 
the  extremal  5:10  ratio  where  the  error  approaches  as  much 
as  25%.  Delay  time  is  more  sensitive  to  k  variations  than 
to  changes  in  B  but  again  errors  exceed  0.03  seconds  only  at 
a  5:10  ratio  and  only  reach  .06  at  that  extreme.  The  larger 
k  values  have  much  reduced  biasing  effect  because  of  the 
filtering  effects  of  the  instrument  response. 

In  total,  the  effect  of  assuming  the  same  source  time 
function  seems  to  not  produce  significant  bias  to  our  es¬ 
timates  of  pP  amplitudes  and  delays.  In  light  of  those 
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Figure  12.  Inversion  solutions  for  synthetic  problem  investigating 
the  effect  of  changes  in  rise  time  parameter  K  for 
Von  Seggern-Blanford  source.  Reference  event  has  K=20, 
B=1  but  K  varying  from  5  to  15.  Only  for  extreme 
model  (K=5)  does  significant  bias  occur. 
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Figure  13: 

Inversion  solutions  for  synthetic  problem  investigating  the 

effects  of  changes  in  overshoot  parameter  B  for  Von  Seggern-  | 

Blandford  source.  Reference  event  has  K=10,  B=l.  Other  event 
has  K=10  but  B  varying  from  0  to  3. 


t 


47 


findings,  and  of  the  results  of  the  other  parameter  studies 

9  discussed  earlier,  we  can  approach  the  application  of  the 

technique  to  observed  data  with  confidence.  The  results  of 
the  application  of  the  method  to  CHTO  and  ANMO  recordings  of 

»  several  Eastern  Kazakh  events  are  discussed  in  the  following 

section. 

t 

t 
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IV.  APPLICATION  TO  DATA 

A  suite  of  seven  presumed  underground  nuclear  ex¬ 
plosions  in  Eastern  Kazakh,  recorded  at  SRO  stations  in 
Thailand  (CHTO)  and  New  Mexico  (ANMO),  were  chosen  as  the 
data  base  for  the  initial  application  of  the  relative 
waveform  technique.  Good  quality  digital  recordings  for  all 
seven  were  events  available  at  both  observatories.  In 
addition,  the  choice  of  these  two  stations  provides  a  range 
in  distance  (~35°  to  CHTO,  ^95°  to  ANMO),  in  azimuth  (south¬ 
east  to  CHTO,  north  to  ANMO)  from  the  source  region,  and  in 
receiver  function. 

The  seismograms  at  the  two  stations  for  all  seven 
events  are  shown  in  Figure^  14  -and— ±5,  with  all  traces 
normalized  to  unit  maximum  amplitude.  These  figures  illus¬ 
trate  the  fact  that  the  effects  in  short  period  seismograms 
of  differing  attenuation  and  near  receiver  earth  structure 
at  two  stations  can  be  considerably  larger  than  the  effects 
of  differing  source  time  functions  and  burial  depths  for  two 
events  recorded  at  a  single  station.  The  ANMO  seismograms 
clearly  show  greater  complexity  than  do  the  CHTO  seis¬ 
mograms,  a  result  of  a  more  complicated  receiver  function. 
In  addition,  the  ANMO  seismograms  clearly  show  greater 
attenuation  than  do  the  CHTO  seismograms. 

The  seven  seismograms  recorded  at  CHTO  show  relatively 
subtle  differences.  These  differences  are  primarily  the 
ratio  of  first  trough  to  second  peak  amplitude,  width  of  the 
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Figure  14: 

Data  seismograms  for  all  events  used  in  this  study  at  stations 
CHTO  and  ANMO.  All  traces  scaled  to  unit  maximum  amplitude. 
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second  peak,  and  amplitude  of  the  second  trough.  These 
differences  are  consistent  in  that  small  second  trough 
amplitude  correlates  with  a  narrow  second  peak.  This  is  the 
type  of  effect  we  expect  to  see  if  the  waveform  differences 
are  caused  by  changes  in  pP  time  and  amplitude.  It  is  also 
worthy  of  note,  at  this  point,  that  two  pairs  of  seismograms 
at  CHTO  show  striking  similarity.  Events  5  and  8  have 
virtually  identical  seismograms,  as  do  events  6  and  7. 
These  event  pairs  are  shown  in  Figures  15  and  16. 

The  greater  complication  of  the  AMMO  receiver  function 
unfortunately  prevents  identification  of  simple,  correlated 
differences  between  events  as  was  apparent  at  CHTO.  Fur¬ 
ther,  it  is  of  some  interest  to  note  that  greater  differ¬ 
ences  exist  for  event  pairs  5  and  8,  and  6  and  7,  at  ANMO 
than  were  apparent  at  CHTO. 

One  method  of  applying  the  relative  waveform  inversion 
technique  is  a  simultaneous  inversion  of  all  seven  events. 
An  alternative  method  is  to  use  a  weighted  regression  tech¬ 
nique  on  the  results  of  two  and  three  event  inversions,  with 
the  weightings  determined  by  the  model  variance  matricies 
obtained  for  each  inversion  run.  The  latter  method  was  used 
in  this  study. 

In  evaluating  the  behavior  of  the  relative  waveform 
inversion  technique,  it  was  found  that  the  absolute  pP 
amplitudes  and  delay  times  were  generally  less  well- 
constrained  than  the  differences  in  amplitudes  and  delay 
times  between  different  events.  This  suggests  that  more 
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Figure  15: 

Event  seismograms  5  and  8  are  shown  above  for  station 
CHTO.  These  two  events  from  Eastern  Kazakh  have  almost 
identical  waveshape  even  out  to  6  seconds.  They  are 
however,  very  different  from  events  6  and  7  in  both 
pulse  width  and  second  trough  depth. 
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Figure  16.  Event  seismograms  6  and  7  are  shown  above  for  station 
CHTO.  Both  events  occurred  in  Eastern  Kazakh  and 
appear  to  be  almost  identical  at  this  S.E.  Azimuth. 
The  seismograms  have  a  narrow  width  and  a  shallow 
second  trough  which  was  uncharacteristic  of  our  other 
data. 
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Figure  17.  Event  seismograms  6  and  7  are  shown  above  for  station 
ANMO.  Both  events  occurred  in  Eastern  Kazakh.  Unlike 
the  CHTO  seismograms  for  the  same  events,  these  wave¬ 
forms  are  noticably  different  in  relative  peak  heights 
shape.  This  may  be  an  indication  of  an  azimuthally 
dependent  source. 
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stable  estimates  of  amplitude  and  delay  time  could  be  ach¬ 
ieved  by  re-parameterizing  the  inversion  scheme  to  solve  for 
the  differences  and  sums  of  those  parameters  and  then  emp¬ 
loying  a  weighted  regression  to  determine  best  estimates  of 
the  parameters  themselves .  This  is  the  approach  we  have 
adopted  in  our  present  analysis. 

Having  once  inverted  the  observed  seismograms  for  the 
difference  parameters,  the  best  estimates  of  the  absolute 
amplitudes  and  delays  were  obtained  by  solving  the  system  of 
equations  expressed  by: 


where 


£  c  x  =  £  b 

b  =  data  column  vector  = 


I 

i 


x  =  best  estimate  vector  = 


c  =  a  basis  operator  matrix  which  maps  x 
into  b,  and 

£  =  a  weighting  vector 

is  the  parameter  difference  between  events  i  and  j, 
Z^j  is  the  parameter  sum  between  events  i  and  j;  A^ 
and  are  the  best  estimates  of  A^  and  respectively. 
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The  weighting  factors  are  proportional  to  the  reciprocal  of 
the  variance.  As  a  consequence,  the  differences  )  are 
weighted  more  heavily  than  the  sums  (I  —  ). 

The  errors  from  this  approach  may  be  determined  by  the 
weighted  standard  deviations.  Fifty  percent  error  bounds 
are  then  given  by: 

*P(A-A  )2 

eA  =  .67  (N-l )Zp  ,  and 


e  =  .67 
T 

Reduction  of  relative  direct  P  data  (S  =  c-/c- )  for 

*  J 

each  inversion  sum  is  done  in  the  same  manner  except  S  is 
first  linearized  using  the  equation: 

Si;j  =  In  (c^Cj)  =  In  ci  -  In  c^ 

Best  estimates  of  In  c^,  expressed  as  In  ci  are  obtained  by 
using  the  weighting  scheme  discussed  above.  Since  S—  is  a 
relative  measure,  we  need  a  relative  error  estimate.  The 
error  in  for  event  i  relative  to  event  j  may  be  ex¬ 

pressed  as: 


Ip(T-T) 

(N-l)Ip 
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£  =±(sij  -§ij>  '  sij  (l-«±Rij) 

Where  is  the  residual  defined  by  the  matrix  equation: 

R  =  b  -  C  x 

This  error  estimate  relative  to  event  j  includes  ran¬ 
dom  and  small  systematic  components.  Large  scale  mis- 
modeling  due  to  bias  or  extra  arrivals  cannot  be  resolved  in 
this  manner. 

In  our  experiments,  two  sets  of  results  were  obtained 
from  the  available  data.  In  the  first  stage,  only  station 
CHTO  was  employed.  Difficulties  with  parameter  resolution 
were  encountered  for  some  data  pairs,  however,  which  compli¬ 
cated  the  interpretation  of  the  results.  The  problem  was 
corrected  by  adding  additional  stations  to  the  inversion. 
Better  resolution  could  also  be  obtained  by  adding  another 
seismogram  to  the  data  base.  The  second  set  of  results  were 
obtained  from  2  or  3  station  inversions  and  are  considered 
to  be  more  reliable. 

Table  1  lists  the  second  set  of  reduced  pP  amplitudes 
and  arrival  times  for  each  event,  with  50%  error  bounds  as 
derived  from  the  error  estimation  methodology  outlined 
above.  Figure  18  is  a  graph  of  the  results  of  the  second 
round  of  calculations  using  both  stations  plotting  amplitude 
vs.  arrival  time  (and  including  calculated  error  bounds). 
Also  plotted  are  the  single  station  results  (indicated  by 

SGI-R-81-048 


58 


Figure  18.  Inversion  results,  giving  effective  pP  amplitudes 
and  delay  times  relative  to  direct  arrival.  Events 
divide  into  two  populations  by  pP  amplitude.  Rel¬ 
ative  solutions  for  large  pP  events  show  little 
change  between  one  and  two  station  inversions,  while 
small  pP  events  show  larger  changes.  Seismograms 
shown  are  comparison  of  data  with  synthetics  for 
CHTO .  Synthetics  contain  inversion  solutions  for 
pP  time  and  amplitude,  von  Seggern-Blanford  time 
function  with  k=10  B=1  and  for  CHTO. 
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small  circles).  No  error  bounds  were  computed  for  the 
t  single  station-derived  model  parameters. 

A  comparison  of  these  two  sets  of  results  illustrates 
several  notable  points.  In  both  sets  of  results  two  popula- 
»  tions  of  events  are  evident.  The  first  population,  con¬ 

sisting  of  events  3,  5,  8  and  9,  have  a  large  (~.7-.9)  pP 
arrival.  Differences  between  pP  amplitudes  and  arrival 
9  times  among  these  events  are  nearly  identical  for  the  one 

station  and  two  station  inversions.  There  is  some  change  in 
absolute  pP  amplitude  between  the  one  and  two  station  re- 
t  suits.  This  is  probably  just  a  reflection  of  the  fact  that 

absolute  pP  amplitudes  are  much  less  well  constrained  by  the 
relative  waveform  inversion  than  are  relative  pP  amplitudes. 
4  Some  possibility  does  exist,  however,  that  the  overall  shift 

in  pP  amplitude  is  due  to  systematic  azimuthal  variation 
either  in  pP  amplitude  or  in  amplitude  of  an  additional, 

«  unmodeled  arrival . 

The  second  population,  consisting  of  events  2,  6,  and  7 

has  pP  amplitudes  considerably  smaller  than  expected  from 
I  linear  elastic  theory.  We  note  that  both  pP  amplitudes  and 

times  change  for  tese  events  between  the  one  and  two  station 
inversions,  not  only  in  an  absolute  sense  but  relative  to 
9  other  events  as  well.  It  is  to  be  expected,  perhaps,  that 

arrival  times  will  be  poorly  constrained  for  small  arrivals 
(for  zero  amplitude  the  arrival  time  is  totally  uncon- 
t  strained),  and  this  result  in  in  fact  observed  in  the  model 

variance  matrices  for  the  small  pP  population.  The  changes 
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in  pP  amplitud  between  the  one  and  two  event  inversions, 
however,  cannot  be  so  easily  explained  away.  Nor  can  the 
fact  that  events  6  and  7  are  virtually  identical  at  CHTO, 
but  show  more  substantial  differences  at  ANMO.  This  sug¬ 
gests  either  the  presence  of  a  sec.  •  i.  probably  positive, 
arrival  with  different  azimuthal  behavior  for  each  of  these 
events,  or  that  whatever  effect  suppresses  pP  for  these 
events  shows  substantial  azimuthal  variation.  The  con¬ 
sistency  of  relative  times  and  amplitudes  in  the  first 
population  leads  us  to  favor  the  second  explanation,  al¬ 
though  considerably  more  work  is  clearly  required  in  this 
area. 

By  design,  the  relative  waveform  inversion  technique 
explicitly  accounts  for  the  interference  of  secondary  arri¬ 
vals  when  determining  overall  signal  amplitude.  Hence,  such 
amplitude  determinations  may  be  more  consistent  than  such 
standard  measures  of  signal  level  as  mfc  or  the  amplitude  of 
selected  half-cycles  of  the  short  period  waveform.  In  order 
to  evaluate  this  effect,  the  relative  amplitudes  of  the 
direct  phases,  as  determined  by  the  relative  waveform  in¬ 
version,  have  been  compared  to  measurements  of  the  A-phase 
(first  peak),  B-phase  (first  peak  to  first  trough)  and 
C-phase  (first  trough  to  second  peak)  measured  amplitudes  of 
all  events  with  respect  to  event  #7.  Figures  19  and  20  plot 
the  results  of  these  comparisons  for  stations  CHTO  and  ANMO 
respectively.  A  fairly  good  consistency  is  achieved  at 
station  CHTO  but  the  near-receiver  complexity  at  ANMO  in¬ 
duces  large  scatter  in  the  measured  amplitudes.  Only  the 
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Figure  19.  Comparison  of  inversion  amplitude  with  other 
amplitude  measures  for  all  events  at  station 
CHTO .  Simplicity  of  receiver  function  and  good 
signal  to  noise  ratio  give  reasonably  consistent 
results  for  all  measures  used.  All  measurements 
referenced  to  event  #7. 


STATION  CHTO  AMPLITUDE  RATIOS  TAKEN  FROM 


to 

t~> 

z 

u 

X 

w 

CS 

3 

CO 

< 

w 


E-i 

U 

w 

OS 


Q 

2 

< 

to 

s 

3 

to 

w 

OS 


US 

< 

w 

2 

E-* 

to 

OS 

M 

b 

Q 

2 

< 

W 

to 

< 

2 

2 

I 

u 


w 

to 

< 

2 

2 


03 


b 

O 


I - - - 1  <C 


63 


us 

< 

z 

w 

o 

2 

M 

« 

M  CO 

to 

co 

2  CH 

Eh 

< 

< 

2  3 

CO 

2 

2 

U  3 

2 

2 

2 

>  CO 

i— i 

1 

1 

z  w 

b 

04 

u 

M  2 

1 

< 

1 

0 

1 

O 

1 

h-+- 

2 

W 

> 

w 


w 

u 

z 

w 

2 

W 

b 

tJ 

2 


oc 

1  at 


r- 

=«= 


h®o0 


2 

fcj 

VO  CO 

-=tt  z 

3 
Z 

E- 

z 

w 

> 

« 

in 

'  =it 


Mo  *0 


NS  o 


0) 

u 

3 

O' 

•H 

b 


Figure  20.  Comparison  of  inversion  amplitudes  with  other 

amplitudes  measures  for  ANMO.  Inversion  results 
are  most  consistent  with  C-amplitudes .  Various 
measures  give  different  relative  amplitude  due 
to  receiver  function  complications  and  noise. 
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B-phase  measurements  at  ANMO  show  something  approaching  a 
stable  estimate.  Finally,  a  comparison  was  made  between  the 
final  inversion  amplitudes  ratios  and  the  m^  amplitude 
ratios  as  reported  in  the  Preliminary  Determination  of 
Epicenters  complied  by  the  U.S.G.S.  We  note  at  this  time 
that  PDE  magnitude  determinations  often  show  errors  of  a  few 
tenths  of  m^  units.  Results  of  this  comparison  are  shown  in 
Figure  21,  with  both  sets  of  amplitudes  referenced  to  event 
#7.  It  is  interesting  to  note  that  for  events  with  large  pP 
arrivals,  mb  based  amplitude  determinations  show  these 
events  to  be  much  larger  compared  to  event  #7  than  do  the 
inversion  results.  This  effect  is  not  observed  for  the 
small  pP  population.  A  comparison  of  the  inversion  results 
with  a  more  carefully  done  m^  study,  presented  in  Part  II  of 
this  report,  shows  the  effect  much  more  strikingly.  Results 
of  this  study  indicate  that  the  bias  in  relative  yield 
estimation  between  events  having  small  and  large  pP  arrivals 
can  reach  40-50%,  with  the  amplitude  of  the  event  with  large 
pP  being  overestimated  by  m^  relative  to  the  event  with 
small  pP. 
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Figure  21.  Comparison  of  source  amplitude  ratios  as 

determined  by  inversion  and  by  M^.  Magnitude 
are  from  P.D.E.  All  events  are  referenced  to 
event  #7. 
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V.  CONCLUSIONS 

The  application  of  the  relative  waveform  inversion 
technique  thus  far  has  clearly  demonstrated  several  im¬ 
portant  features.  First,  it  is  possible  to  obtain  well- 
resolved,  stable  estimates  of  relative  event  size,  pP  am¬ 
plitude,  and  pP  delay  times  from  relatively  sparse  data  sets 
and  if  necessary,  from  noisy  or  complex  observations. 
Moreover,  the  biasing  effects  of  fairly  small  changes  in  pP 
amplitude  and  delay  time  have  been  demonstrated  on  con¬ 
ventional  yield  estimation  methodologies.  This  biasing 
seems  generally  to  result  in  over-estimates  (from  m^-type 
measures)  of  the  yield  of  underground  explosions  of  events 
with  large  pP. 

This  technique  shows  great  promise  in  yield  estimation 
and  depth  estimates  of  buried  explosions,  but  clearly  sub¬ 
stantially  more  research  efforts  are  required.  The  data 
base  examined  to  date  represents  a  relatively  small  number 
of  events,  azimuths  from  source  to  receiver,  and  stations. 

A  more  exhaustive  testing  program  should  be  undertaken  in 
order  to  fully  test  the  procedure  as  well  as  provide  a  basis 
for  any  decision  to  incorporate  this  kind  of  analysis  into 
formal  monitoring  procedures.  Further  research  is  alo 
necessary  on  the  incorporation  of  relative  waveform  methods 
with  absolute  waveform  methods  to  obtain  even  more  accurate 
absolute  amplitude  and  delay  time  measurements. 
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APPENDIX 

Partial  derivatives  for  two  seismogram  relative  waveform  inversion 

The  source  model  for  the  two  seismogram  case  is  given 
by 

S,  =  1/c  (6(t)  -  a.  6 > t-t,  )) 

1  11/  (A.i) 

S2  =  c  (6(t+t0)  -  a2  6  ( t+t 0— x 2 )) 

The  error  function  is  given  by 

e  =  (fl  *  s2  -  S1  *  f2>  *  W 

=  *  s2  “  Si  *  F2  where  (A. 2) 

Fi  =  fi  *  W' 

Then,  defining  m j ,  the  model  vector,  by 
m  —  ( a^ ,  a2 ,  c,  x  ^ ,  ^i*  ^2^ 

Letting  e.  be  the  jth  time  point  in  the  time  series  e 
we  find 

A- •  =  9ei  is  given  by 

3  3m  j 

3e 

i  =  IF  (t  -  t  ) 

3a1  c  2  i  1 

■  cFl(ti  +  -  l2> 

3a„ 


